function [f,x,xp,y,yp,symparams,Phi] = DSGEmodel_yieldCurve(unitFree)
% This function reports the equations for the New Keynesian model with Calvo prices
% in f, and it reports the state vector (x and xp) and the control vector (y and yp)

%Define parameters
syms BETTA B CHI CHI0 THETA DELTA ALFA PHI PHIzero ZI KAPAw ETA RHOR PHIpai PHIpai_1 ...
     PHIy PHIy_1 PHIc PHIc_1 PHIl NU U0 U0d...
     OMEGAz OMEGAd OMEGAn OMEGAp OMEGAa ...
     RHOz  RHOd   RHOn  RHOp RHOa ...
     Kss OUTPUTss PAIss MUZss Css AA lss Wss;
Rss = sym('Rss');
symparams=[BETTA,B,CHI,CHI0,THETA,DELTA,ALFA,PHI,PHIzero,ZI,KAPAw,ETA,...
     RHOR,PHIpai,PHIpai_1,PHIy,PHIy_1,PHIc,PHIc_1,PHIl,NU,U0,U0d,...
     OMEGAz,OMEGAd,OMEGAn,OMEGAp,OMEGAa,...
     RHOz ,RHOd   RHOn  ,RHOp,RHOa, ...
     Kss OUTPUTss PAIss MUZss Css AA lss,Rss,Wss];

%Define variables   
syms c_cu   r_cu    pai_cu   w_cu  l_cu  output_cu  mc_cu  evf_cu  vf_cu...
     c_cup  r_cup   pai_cup  w_cup l_cup output_cup mc_cup evf_cup vf_cup...
     c_ba1  muz_cu   d_cu  n_cu  paistar_cu  a_cu...
     c_ba1p muz_cup  d_cup n_cup paistar_cup a_cup;

% Bond prices
syms       p1_cu   p2_cu   p3_cu   p4_cu   p5_cu   p6_cu   p7_cu   p8_cu   p9_cu   p10_cu ...
           p11_cu  p12_cu  p13_cu  p14_cu  p15_cu  p16_cu  p17_cu  p18_cu  p19_cu  p20_cu ...
           p21_cu  p22_cu  p23_cu  p24_cu  p25_cu  p26_cu  p27_cu  p28_cu  p29_cu  p30_cu ...
           p31_cu  p32_cu  p33_cu  p34_cu  p35_cu  p36_cu  p37_cu  p38_cu  p39_cu  p40_cu;
syms       p1_cup  p2_cup  p3_cup  p4_cup  p5_cup  p6_cup  p7_cup  p8_cup  p9_cup  p10_cup ...
           p11_cup p12_cup p13_cup p14_cup p15_cup p16_cup p17_cup p18_cup p19_cup p20_cup ...
           p21_cup p22_cup p23_cup p24_cup p25_cup p26_cup p27_cup p28_cup p29_cup p30_cup ...
           p31_cup p32_cup p33_cup p34_cup p35_cup p36_cup p37_cup p38_cup p39_cup p40_cup ...    

%% The equations in f
% EQ1: The expression for the value function (minus)
%vf_cup = - (U0 + d_cup*(1/(1-CHI)*((c_cup-B*c_cu/muz_cup)/Css^CHI0)^(1-CHI)+U0d + n_cup*PHIzero/(1-1/PHI)*(1-l_cup)^(1-1/PHI)) ) + BETTA*AA*evf_cup^(1/(1-ALFA));    
if unitFree == 1
    f1 =  -1 + (-(U0 + d_cu*(1/(1-CHI)*((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(1-CHI)+U0d + n_cu*PHIzero/(1-1/PHI)*(1-l_cu )^(1-1/PHI)) ) ...
                + BETTA*AA*evf_cu ^(1/(1-ALFA)))/vf_cu;
else
    f1 =  -vf_cu - (U0 + d_cu*(1/(1-CHI)*((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(1-CHI)+U0d + n_cu*PHIzero/(1-1/PHI)*(1-l_cu )^(1-1/PHI)) )  + BETTA*AA*evf_cu ^(1/(1-ALFA));
end

% EQ 2: Household's First-order condition for labour
if unitFree == 1
    f2 = -1 + (KAPAw*Wss+(1-KAPAw)*(n_cu*PHIzero*(1-l_cu)^(-1/PHI)/((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(-CHI)*Css^CHI0))/w_cu;
else
    f2 = -w_cu + KAPAw*Wss+(1-KAPAw)*(n_cu*PHIzero*(1-l_cu)^(-1/PHI)/((c_cu -B*c_ba1/muz_cu)/Css^CHI0)^(-CHI)*Css^CHI0);
end

% EQ 3: Euler-equation for the one-period interest rate
% The real stochastic discount factor
mReal_cup = BETTA*d_cup/d_cu*...
    (AA*evf_cu^(1/(1-ALFA))/(vf_cup*muz_cup^((1-CHI)*(1-CHI0)) ))^ALFA...   
    *(c_cup-B*c_cu/muz_cup)^(-CHI)/((c_cu-B*c_ba1/muz_cu)^(-CHI))*muz_cup^(-CHI*(1-CHI0)-CHI0);  
if unitFree == 1
    f3 = -1 + mReal_cup/pai_cup*r_cu;
else
    f3 = -1 + mReal_cup/pai_cup*r_cu;
end

% EQ 4: Firms FOC for labour
if unitFree == 1
    f4 = -1  + (w_cu/((1-THETA)*a_cu *Kss^THETA*l_cu ^-THETA))/mc_cu;
else
    f4 = -mc_cu  + w_cu/((1-THETA)*a_cu *Kss^THETA*l_cu ^-THETA);
end

% EQ 5: First-order condition for prices
if unitFree == 1
    f5 = ((1-ETA)+ZI*mReal_cup*(pai_cup/PAIss^NU-1)*pai_cup/PAIss^NU*output_cup/output_cu*muz_cup ...
        - ZI*(pai_cu/PAIss^NU-1)*pai_cu/PAIss^NU)/(ETA*mc_cu) + 1;
else
    f5 = (1-ETA)+ZI*mReal_cup*(pai_cup/PAIss^NU-1)*pai_cup/PAIss^NU*output_cup/output_cu*muz_cup ...
        + ETA*mc_cu - ZI*(pai_cu/PAIss^NU-1)*pai_cu/PAIss^NU;
end

% Taylor rule
f6= -r_cu + Rss*exp((PHIpai*(log(pai_cu)-(log(PAIss) + log(paistar_cu)) )...
         + PHIpai_1*(log(pai_cup)-(log(PAIss) + log(paistar_cup))) ...
         + PHIy*log(output_cu/OUTPUTss)...
         + PHIy_1*log(output_cup/OUTPUTss)...
         + PHIc*(c_cu-c_ba1+log(muz_cu)-log(MUZss)) ...
         + PHIc_1*(c_cup-c_ba1p+log(muz_cup)-log(MUZss)) ...
         + PHIl*(log(l_cu/lss))));

    
%EQ7: Agg. ressource constraint
if unitFree == 1
    f7 = - output_cu*(1-ZI/2*(pai_cu /PAIss^NU-1)^2)/(c_cu  + DELTA*Kss) + 1;
else
    f7 = - output_cu*(1-ZI/2*(pai_cu /PAIss^NU-1)^2) + (c_cu  + DELTA*Kss);
end

% EQ 8: The aggregated production function to get output. 
if unitFree == 1
   f8 = -1 + a_cu*Kss^THETA*l_cu ^(1-THETA)/output_cu;   
else
   f8 = -output_cu + a_cu *Kss^THETA*l_cu ^(1-THETA);  
end

% The expected value of the value function 
if unitFree == 1
    f9 = -1 + (vf_cup*muz_cup^((1-CHI)*(1-CHI0))/AA)^(1-ALFA)/evf_cu;
else
    f9 = -evf_cu + (vf_cup*muz_cup^((1-CHI)*(1-CHI0))/AA)^(1-ALFA);
end

%% Nominal Bond Prices 
f10= -1 + (1/r_cu)/p1_cu;
f11= -1 + mReal_cup/pai_cup*p1_cup/p2_cu;
f12= -1 + mReal_cup/pai_cup*p2_cup/p3_cu;
f13= -1 + mReal_cup/pai_cup*p3_cup/p4_cu;
f14= -1 + mReal_cup/pai_cup*p4_cup/p5_cu;
f15= -1 + mReal_cup/pai_cup*p5_cup/p6_cu;
f16= -1 + mReal_cup/pai_cup*p6_cup/p7_cu;
f17= -1 + mReal_cup/pai_cup*p7_cup/p8_cu;
f18= -1 + mReal_cup/pai_cup*p8_cup/p9_cu;
f19= -1 + mReal_cup/pai_cup*p9_cup/p10_cu;
f20= -1 + mReal_cup/pai_cup*p10_cup/p11_cu;
f21= -1 + mReal_cup/pai_cup*p11_cup/p12_cu;
f22= -1 + mReal_cup/pai_cup*p12_cup/p13_cu;
f23= -1 + mReal_cup/pai_cup*p13_cup/p14_cu;
f24= -1 + mReal_cup/pai_cup*p14_cup/p15_cu;
f25= -1 + mReal_cup/pai_cup*p15_cup/p16_cu;
f26= -1 + mReal_cup/pai_cup*p16_cup/p17_cu;
f27= -1 + mReal_cup/pai_cup*p17_cup/p18_cu;
f28= -1 + mReal_cup/pai_cup*p18_cup/p19_cu;
f29= -1 + mReal_cup/pai_cup*p19_cup/p20_cu;
f30= -1 + mReal_cup/pai_cup*p20_cup/p21_cu;
f31= -1 + mReal_cup/pai_cup*p21_cup/p22_cu;
f32= -1 + mReal_cup/pai_cup*p22_cup/p23_cu;
f33= -1 + mReal_cup/pai_cup*p23_cup/p24_cu;
f34= -1 + mReal_cup/pai_cup*p24_cup/p25_cu;
f35= -1 + mReal_cup/pai_cup*p25_cup/p26_cu;
f36= -1 + mReal_cup/pai_cup*p26_cup/p27_cu;
f37= -1 + mReal_cup/pai_cup*p27_cup/p28_cu;
f38= -1 + mReal_cup/pai_cup*p28_cup/p29_cu;
f39= -1 + mReal_cup/pai_cup*p29_cup/p30_cu;
f40= -1 + mReal_cup/pai_cup*p30_cup/p31_cu;
f41= -1 + mReal_cup/pai_cup*p31_cup/p32_cu;
f42= -1 + mReal_cup/pai_cup*p32_cup/p33_cu;
f43= -1 + mReal_cup/pai_cup*p33_cup/p34_cu;
f44= -1 + mReal_cup/pai_cup*p34_cup/p35_cu;
f45= -1 + mReal_cup/pai_cup*p35_cup/p36_cu;
f46= -1 + mReal_cup/pai_cup*p36_cup/p37_cu;
f47= -1 + mReal_cup/pai_cup*p37_cup/p38_cu;
f48= -1 + mReal_cup/pai_cup*p38_cup/p39_cu;
f49= -1 + mReal_cup/pai_cup*p39_cup/p40_cu;

%% The links
% For lagged consumption
if unitFree == 1
    f50 = -1 + c_cu/c_ba1p;
else
    f50 = -c_ba1p + c_cu;
end

%% The linear shocks
%For non-stationary technology
f51 = -log(muz_cup/MUZss) + RHOz*log(muz_cu/MUZss);

%The preference shock
f52 = -log(d_cup) + RHOd*log(d_cu);

%For labor shocks
f53 = -log(n_cup) + RHOn*log(n_cu);

%Inflation target
f54 = -log(paistar_cup) + RHOp*log(paistar_cu);

%stationary technology shocks
f55 = -log(a_cup) + RHOa*log(a_cu);


%Creating function f
f =[ f1; f2; f3; f4; f5; f6; f7; f8; f9;f10;...
    f11;f12;f13;f14;f15;f16;f17;f18;f19;f20;...
    f21;f22;f23;f24;f25;f26;f27;f28;f29;f30;...
    f31;f32;f33;f34;f35;f36;f37;f38;f39;f40;...
    f41;f42;f43;f44;f45;f46;f47;f48;f49;f50;...
    f51;f52;f53;f54;f55];

% Define the variables which needs to be log-transformed
y  = [ c_cu   r_cu    pai_cu  w_cu    l_cu    output_cu  mc_cu  evf_cu  vf_cu...
      p1_cu   p2_cu   p3_cu   p4_cu   p5_cu   p6_cu   p7_cu   p8_cu   p9_cu   p10_cu ...
      p11_cu  p12_cu  p13_cu  p14_cu  p15_cu  p16_cu  p17_cu  p18_cu  p19_cu  p20_cu ...
      p21_cu  p22_cu  p23_cu  p24_cu  p25_cu  p26_cu  p27_cu  p28_cu  p29_cu  p30_cu ...
      p31_cu  p32_cu  p33_cu  p34_cu  p35_cu  p36_cu  p37_cu  p38_cu  p39_cu  p40_cu ];
yp = [c_cup   r_cup   pai_cup w_cup   l_cup   output_cup  mc_cup  evf_cup  vf_cup...
      p1_cup  p2_cup  p3_cup  p4_cup  p5_cup  p6_cup  p7_cup  p8_cup  p9_cup  p10_cup ...
      p11_cup p12_cup p13_cup p14_cup p15_cup p16_cup p17_cup p18_cup p19_cup p20_cup ...
      p21_cup p22_cup p23_cup p24_cup p25_cup p26_cup p27_cup p28_cup p29_cup p30_cup ...
      p31_cup p32_cup p33_cup p34_cup p35_cup p36_cup p37_cup p38_cup p39_cup p40_cup ];
x  = [c_ba1  muz_cu   d_cu  n_cu  paistar_cu  a_cu];
xp = [c_ba1p muz_cup  d_cup n_cup paistar_cup a_cup];


% Make f a function of the logarithm of the state and control vector
f = subs(f, [x,y,xp,yp], exp([x,y,xp,yp]));

Phi = [];
end